The critical exponents of the two-dimensional Ising spin glass revisited: 
Exact Ground State Calculations and Monte Carlo Simulations 
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The critical exponents for T ^ of the two-dimensional Ising spin glass model with Gaussian 
couplings are determined with the help of exact ground states for system sizes up to L = 50 and 
by a Monte Carlo study of a pseudo-ferromagnetic order parameter. We obtain: for the stiffness 
exponent y{= 9) — —0.281 ± 0.002, for the magnetic exponent S = 1.48 ± 0.01 and for the chaos 
exponent ( — 1.05 ±0.05. From Monte Carlo simulations we get the thermal exponent i/ = 3.6 ±0.2. 
The scaling prediction y = —l/i/ is fulfilled within the error bars, whereas there is a disagreement 
with the relation y — 1 — 5. 

PACS numbers: 75.40, 05.45, 75.10 



I. INTRODUCTION 

It is now widely believed that the bond-disordered two- 
dimensional Ising spin glass model with short range in- 
teractions does not hajue. a phase transition at any non- 
vanishing temperatureoo. At zero temperature the spin 
glass is in its ground state (i.e. the spin configuration with 
the lowest possible energy), which might be degenerate 
or unique depending on the probability distribution of 
the spin interactions. This ground state is unstable with 
respect to thermal fluctuations and any non-vanishing 
temperature destroys this long range spin glass order. By 
decreasing the temperatures on the other hand the spa- 
tial correlations grow resulting in a divergence of the spin 
glass susceptibility at zero temperature. This scenario is 
characterized by a set of critical exponents that depend 
on certain features of the bond distribution. Experiments 
on Rb2Cui_2;Co2;F4 clearly confirmed this picturea and 
reported values for the critical exponents, which are com- 
patible with those predicted by the numerical investiga- 
tions. 

The latter has been pursued in four different waira: 
Monte Carlo simulations at finiljfi temperaturesu'El, 
high temperature series expansionQ, transfer matrix 
calculationa3 eI and exact determinaiinn- of ground 
states via cptnbinatorial optimizationt3~Lj or replipa 
optimizationll3. A scaling theory by Bray and MooretZl 
establishes relations between exponents quantifying the 
stiffness of the ground state and the critical exponents 
characterizing the temperature dependent divergence of 
various thermodynamic quantities like correlation length 
or susceptibility. 

With the most recent numerical studies a controversy 
arisen on the critical exponents of the two-dimensional 
Ising spin glass with Gaussian couplings: A Monte Carlo 



study by Liangi (using a Swendsen- Wang-type cluster 
algorithm) and a numerical transfer matrix calculation 
by Kawashima et alB yield a value for the thermal ex- 
ponent i/p^whirh is significantly different from the early 
estimatesafl^El. Moreover, Kawashima et al£3 also study 
the ground state magnetization of this model in an ex- 
ternal field and report a value for the n^agnetic field ex- 
ponent, which is, using a scaling relatiorJlIl, incompatible 
with the stiffness exponent found in domain wall renor- 
malization group studies. 

This observation and the progress in algorithmic de- 
velopments motivated us to revisit the critical exponents 
of the two-dimensional Ising spin glass model. In this pa- 
per we present a synopsis of a zero-temperature (ground 
state) and a finite-temperature (Monte Carlo) approach 
to estimate the numerical values for these critical expo- 
nents. For the former we report results obtained from ex- 
act ground states for the largest system sizes possible to 
date, resulting in the most reliable estimates for the stiff- 
ness exponent y and magnetic field exponent S reported 
so far. For the Monte-Carlo simulations we propose a 
pseudo-ferromagnetic order parameter that is defined by 
a projection of spin configurations onto the exactly know 
ground state and show that the thermal exponent v is 
identical t&nthe values that have been obtained by Bhatt 
and YoungQ studying the Edwards- Anderson (EA) order 
parameter. In the context of domain growth and non- 
equilibrium dynamics this concept has already been in- 
troduced and proven to be usefuLby direct comparison 
with the so called replica overlapc3. 

The two-dimensional Ising spin glass model with a 
Gaussian distribution of couplings, which we consider 
throughout this paper, is defined by the Hamiltonian 



H — JijSiSj — h Si , Si — ±1 



(1) 



1 



where (ij) denotes all nearest neighbor pairs on a L x L 
square lattice with periodic boundary conditions and the 
random interaction strengths obey a Gaussian probabil- 
ity distribution with mean zero and variance one. The 
parameter h denotes an external magnetic field strength. 

The outline of the paper is as follows: In the next 
section we present our results from exact defect energy 
calculations, which provides us with an estimate for the 
stiffness exponent y. In section III we present a conven- 
tional finite size scaling analysis of Monte Carlo data. In 
contrast to earlier investigations we used exact ground 
state configurations instead of replica systems in order 
to establish an order parameter. Section IV focuses on 



the exact calculation of ground state magnetizations in 
an external field and its finite size scaling properties. Sec- 
tion ^ presents a study of the sensitivity of the ground 
state with respect to slight perturbations of the coupling 
strength. The last section is a summary plus discussion. 



II. DEFECT ENERGY 



A. Scaling Theory 



The scaling theory by Bray and Moore0 starts with 
a coarse grained picture for the spin interactions. It hy- 
pothesizes the following scaling ansatz for an effective 
couphng J{L) among (block) spins on length scale L at 
an infinitesimal temperature: 



j(L) - jiy 



(2) 



where J denotes the variance of the original bond dis- 
tribution. For positive stiffness exponent y (sometimes 
9) the coupling becomes stronger on larger length scales, 
which means that it is harder to flip collectively a con- 
nected set of spins of linear dimension L. Thus thermal 
fluctuations are irrelevant and the spin glass ordered state 
persists at low temperature. A negative exponent y, as 
we expect for d = 2, on the other hand indicates the in- 
stability of the spin glass ground state. In this case the 
spin glass transition takes place only at zero temperature 
and y is related to the thermal exponent v determining 
the divergence of the correlation length ^ ^ T^'^: 

The temperature dependence of the correlation length 
^ near T = can be inferred from equating the two 
energy scales set by the effective coupling constant and 
the temperature. At low temperatures where (|^) holds 
one has then 



and therefore 



y 



-\lv 



(3) 



(4) 



In this way the exponent y, which we calculate in this 
section, has to be compared with v determined in finite 
temperature Monte Carlo simulations discussed in the 
next section. 



B. Algorithm and Results 

The problem of finding a spin configuration with low- 
est energy can be transformed into the problem of find- 
ing a maximum weight cut in a special weighted graph, 
that represents the interaction structure of the spin glass 
systemEiJ. This is known under the name Max-Cut prob- 
lem and is in general a NP-complete problemtil. If the 
graph is planar, as in the two-dimensional case with /ree 
or fixed boundariz. conditions, the problem is solvable in 
polynomial timctJ. If one has periodic or anti-periodic 
boundary conditions or if an external field (representable 
as an extra node to which all other nodes are connected) 
is present the graph is not planar any more even in two 
dimensions. Hence, the situation we are studying here is 
indeed a NP-complete problem. 

The results of this section and of sections |^ and ^ 
are based on the application of a so-called branch & cut 
algorithm to the ground state problemtd. This algorithm 
always finds an exact ground state of the given spin glass 
system. For further details about this algorithm and its 
implementation see ref.E3. An important feature of this 
approach is, that the returned solutions are proved to be 
optimal. Exact ground states of grid sizes up to 100 x 100 
can be determined in a moderate amount of computation 
time. The 100 x 100 instances take between 1.5 and 8 
hours, 4 hours on average. Up to grid sizes of 50 each 
run takes less than 15 minutes. 

The NP-completeness is not a serious problem as long 
as the system sizes L one studies are not in the region 
where the exponential dominates, and for really large L, 
where it does dominate, the exponent is very small. In 
the range L < 32 our empirical CPU-times can be fit- 
ted by a power law (r cx L?'^). For bigger systems there 
enters an exponential term « 1.2^ inside the used size 
range. Note that these are only empirical observations 
but no rigorous bound for the complexity. For compar- 
ison Kawashima and Suzukll3 reported about a replica 
optimization method which approximates ground states 
efficiently. They achieved an average CPU-time t(L) for 
systems of linear size L which can be fitted by a power 
law like t oc iF'-^ inside the same size range (L < 32). 
Although Kawashima and Suzuki only approximate the 
ground states while we always find optimal solutions, 
their CPU-times are similar to ours. On average, their 
biggest systems (32 x 32) took 260 seconds, while we 
needed 160 seconds for 40 x 40 spin glasses (1500 seconds 
for 60 X 60). Kawashima and Suzuki used a VAX 6440, 
our computations were carried out on a SPARC 10/612. 

One can determine the defect energy by investigat- 
ing the sensitij«te-|Of the ground state energy to bound- 
ary conditionglZlllj, which can be quantified by a stiff- 
ness exponent measuring the extra energy of a defect 
line through the whole sample. The block coupling J' is 
then given by J' — ^ (Ep ~ Ea)^, where Ep and Ea are 
the ground state energies of the system under periodic 
and anti-periodic boundary conditions, respectively. We 
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FIG. 1. Defect energy Ed as a function of the system size L 
in a log-log plot. The straight line is a least square fit giving 
the exponent y = —0.281. 

compute this value in the following way. First we solve 
the given spin glass system to optimality under periodic 
boundary conditions, i.e., we find an exact ground state 
configuration ujp and its energy Ep — E{ujp). Then we 
choose two neighboring "columns" of spins and multi- 
ply all couplings by —1 that link these two spin sets. By 
this modification of the couplings we impose anti-periodic 
boundary conditions to the original system. With this 
slightly changed objective function we rerun our branch 
& cut code to find a ground state configuration uja with 
energy Ea= E{uja). 

Due to the very small magnitude of Ea — Ep it is nec- 
essary to have a very large number of samples to obtain 
stable statistics. For each size L < 30 of our L x L spin 
glasses we ran [ ^'^" ] samples. The resulting mean val- 
ues of the defect energies versus the system size L are 
shown in Fig. |l|. A least square fit yields the value 

y = -0.281 ±0.002 ^ = 3.559 ± 0.025 (5) 

The errors are statistical errors only. This estimate 
agrees roughly with less accurate early estimates v = 
2.96 ± 0.22|-|and i' = 4.2 ± 0.5 from transfer matrix 
calculationgLl as well as = 3.56 ± 0.06 and = 3.4 ±0.1 
from domain wall reiiprmalization calculationsla. Note 
that Bray and MooreliZi report an estimate y = —0.291 ± 
0.002, which has an error bar that is identical to ours. 
However, their maximum system size is L = 12 and they 
did not calculate exact ground states. 

Our result for y implies a value for ly, if the scaling 
prediction is corrapiL that differs substantially from 
more recent estimatesQ'EI. Since in these works the ther- 
mal exponent v has been determined directly, we shall 
do this, too, in the next section. 



III. MONTE CARLO RESULTS 

In this section we present our results from finite tem- 
perature Monte Carlo simulations. For this purpose we 
introduce first a number of quantities that are of interest 
for studying the critical properties of Ising spin glasses 
in zero external field h. 



A. Scaling Relations and Methodology 

A characteristic feature of a spin glass transition at a 
temperature Tc (which might be zero) is the divergence 
of the so called spin glass susceptibility 

x = ^ E[(^^^^-)'w, (6) 

<i]> 

where [. . .]av denotes the average over the quenched dis- 
order and (...) a thermal average. Approaching the tran- 
sition temperature from the paramagnetic phase one ob- 
serves X {T — Tc)~^ , which defines the susceptibility 
exponent 7. As already mentioned there is a diverging 
length scale at the transition, the spin glass correlation 
length ^ {T — Tc)~'^, which governs the scaling form of 
the correlation function near Tc'. 

G(r) = [< SA+r >%. ^ r-^^-^+^^-gir/O (7) 

Obviously 7 = (2 — ri)iy. In the case Tc > there 
is a non-vanishing Edwards- Anderson order parameter 
qsA — [{Si)'^]s,v below the transition and one has qsA ~ 
(Tc — T)f^ for T < Tc, the order parameter exponent f3 
obeying the hyperscaling relation f3 = ^{d — 2 + rj). 

In two dimensions Tc = and since we are concerned 
with a continuous bond distribution, in which case the 
ground state is non-degenerate, one has 

?7 = 

/3 = (8) 
7/1. = 2 

Thus we are left with a single unknown exponent j/, which 
we determined from the scaling behavior of the suscepti- 
bility and the Binder cumulant. In contrast to previous 
investigations we used exact ground states to define a 
pseudo-ferromagnetic order parameter via 

M = mu (9) 

with 

1 ^ 

9=]^E^''^° iN = L') (10) 

where Sf denotes the value of the spin at site i in one of 
the two ground state configurations. Note that in con- 
trast to the EA-order parameter here is only one fluctu- 
ating fleld involved, which would in principle reduce the 
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order parameter exponent to /9/2. Since we have Tc = 
and thus /3 = this is not relevant here. The correspond- 
ing order parameter susceptibiUty is defined via 



XL - N[{q^% 



(11) 



For the finite size scaling form of the susceptibility we 
expect according to d§) 



XL 



(T) ^ L^xiL^^'^T). 



(12) 



The second quantity we studied was the disorder aver- 
aged Binder cumulant 



1 

5L- 2 



(13) 



Since this is a dimensionless combination of moments its 
finite size scaling form is 



9l{T) 



(14) 



This quantity provides us with a second independent es- 
timate for the scaling exponent v. 

We applied single spin flip Glauber dynamics to per- 
form our simulations with a spin flip probability given 
by 



w{Si 



-Si) = 



1 



1 + cxp(A£:/r)' 



(15) 



and A_B being the energy difference between the new and 
the old state. Time is measured in Monte Carlo sweeps 
(MCS) through the whole lattice. 
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FIG. 2. Plot of the susceptibility at T = 1.4 for L = 10 
averaged over 192 samples. We compare the values of Xl(^™) 
obtained from systems initialized with a ground state con- 
figuration with those obtained from systems with a random 
initial configuration. Obviously the results become time inde- 
pendent (aside from statistical fluctuations) if both estimates 
agree. 

The estimate of the critical exponents necessitates the 
determination of the equilibrium values of the thermody- 
namic quantities of interest. Due to the slow relaxation of 



spin glasses it is difficult to decide whether the values are 
stationary or not, because it is hard to discriminate be- 
tween real- and quasi-stationary values of the functions. 
This is why we used a definite criterion analogous to the 
criterion introduced by Bhatt and Youngj: 

We simulated two replicas of the system, one which 
has been initialized with a random configuration and the 
other with the ground state configuration. From Fig.|^ 
it can be clearly seen that if both estimates agree we 
obtained a time independent value of the susceptibility, 
which we took as our equilibration criterion. 



B. Results 

We studied the temperature dependent scaling behav- 
ior system sizes extending from L = 6 to L — 12. 
The number of samples is chosen that approximately 
N ■ ^samples — const, holds. We simulated at least 128 
samples for L = 12. 
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FIG. 3. Equilibrium values of the susceptibiUty depending 
on temperature and system size. 

Fig.^ shows the equilibrium values of the susceptibility 
for various system sizes. We could reach the equilibrium 
value of the susceptibility in the chosen time interval for 
T > 1.0 {L > 6) and T > 0.8 {L = 6) respectively. With 
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FIG. 4. Results of gL- Only equilibrium- values are shown. 



this data we got an estimate from the scaling-ansatz (12). 
The best data collapse we obtained (see FigJ|) for 



z/ = 3.4 ±0.2. 



(16) 



The error bars denote the interval of exponents where we 
get an indistinguishable data collapse. 

Fig.^ shows the equilibrium values of gj^ for the same 
samples. The equilibrium value could be reached within 
the same time interval for some smaller temperature. 
This is why the data is more sensitive to changes of the 
value of the critical exponent. Thus we obtained the best 
data collapse for 



V = 3.7±0.1. 



(17) 



in agreement within the error bars with the value de- 
termined above. Concluding we obtain an average value 
of 



3.6 ±0.2 



(18) 



for the critical exponent from our Monte Carlo simu- 
lations. This value agrees well with the estimate (|^) 
that we obtained from the defect energy calculations 
in the last section. It differs substantially from the 
more xecent estimates obtained by a cluster Monte Carlo 
studyQ {v 2.0 ± 0.2) and a numerical transfer matrix 
calculationla (i/ = 2.08 ± 0.01). 



IV. GROUND STATE MAGNETIZATION 



A nonzero external field h induces a non-vanishing 



magnetization m = N ^ 12iLi '^i ^ system 
ground state {5°}. The relat ion between magnetization 
and field strength is highly nontrivial in general and mo- 
tivates the introduction of a new exponent S characteriz- 
ing this relation in the infinite system (L —^ oo) for small 
fields J): 



with 



(19) 



The corresponding finite size scaling form and a scaling 
relation between 6 and the already known exponent y can 
be obtained by the following argumentO: 

If the ground state is non-dcgencrate the spins are ran- 
domly oriented within an infinitesimal field at T = 0. 
Hence the magnetization of a finite system in zero 
field is a random variable with variance l/N, implying 
mL{h = 0) ^ L"'^!'^ . As a further consequence of the ran- 
dom orientations the total magnetic moment of a block 
spin of linear dimension L is of order L"^/^, thus the mag- 
netic field on this length scale has to be rescaled accord- 
ing to 



(20) 



in contrast to a ferromagnet, where we would have 
h(V) ^ L'^h. For nonzero field (at T = 0) one would 
expect mL{h) ■ L'^^^ to be a function of the dimensionless 
ratio of energy scales J{L) and h{L) only, thus 



rriLih) = L-'^/^ m{L'^/^-y hJ-^) 



(21) 



with rh{x — + 0) = const.. Since for L ^ oo the L- 
dependence of the magnetization has to drop out it is 
m{x oo) 2:''/(<i-2K)^ Moreover, in this limit we have 
to recover ([l9| ) , which implies for d — 2 



(22) 



Rewriting ( pT| ) slightly for our our purposes yields 

mL{h) = L-^m{Lh^/^) (23) 

with m{x 0) = const, and rn{x — > oo) oc x. Note 
that (p3|) should hold independently of the correctness 
of the above derivation of the scaling relation (|2|) : The 
length scale induced by the magnetic field is given by 
h~^/^ and (|2^) is simply the finite size scaling form one 
would expect for the magnetization. 

With our branch & cut algorithm we are not only 
able to compute m{S, h) for a sample S for some spe- 
cific values-uf h like other authors did (see Kawashima 
and Suzukilla). We can evaluate the complete piecewise 
constant function 771(5*, h) for each sample. We do this by 
starting at ft, = 0, computing the ground state, and find- 
ing the next increased value of h for which the current 
ground state loses optimality using a sensitivity analysis 
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FIG. 5. Scaling plot for the ground state magnetiza- 
tion: LmhiK) versus U\}^^ for various system sizes with 
1/(5 = 0.675. Note that for high fields h ^ oo the curves 
have to saturate at L ■ mL{h oo) — L. 

technique0. At that point we compute the new ground 
state. We do this up to a given field strength or until 
saturation occurs. 

This technique gives us the (averaged) function mL{h) 
for each system size L with any arbitrary resolution. We 
used systems of sizes L e {10, 20, 30, 40, 50, 60} and com- 
puted the ground states for \^~\ samples for each size L. 
We judged the data collapse in a plot tolL versus Lh^/^ 
as shown in Fig. 5 by visual inspection and using cubic 
spline interpolation. In the figure we have included some 
error bars for the L = 50 and the L = 60 curves to show 
the typical errors. The errors decrease with decreasing 
system size, because of the increasing number of samples. 

We obtained the best data collapse at 



1/(5 = 0.675 ±0.005 



1.481 ±0.011 (24) 



a value that agrees well with the result of the groupd 
state magnetization study by Kawashima and Suzukill3. 
This value together with estimate y — 0.281 (|^) from the 
defect energy calculation implies that the scaling hypoth- 
esis is significantly violated. Since we estimated v 
directly in a Monte Carlo simulation we conclude that it 
also not legitimate to infer from the magnetic exponent 
5 via ( p^ and ( p^ ) that-tie thermal exponent v should 
be close to 2 as found inO'LI. 



V. CHAOS EXPONENT 

One of the peculiar features of spin glasses is their ex- 
treme sensitivity with respect to parameter changesE3, 
like small temperature, field or coupling variations. For 



the ground state properties this means that a slight per- 
turbation of the initial set of couplings leads to a com- 
plete reorganization of the original ground state over a 
length scale that depends on the strength of the per- 
turbatiari_-,This overlap length is expected to behave as 
foUowsHe: 

Let us modify the interactions by replacing each cou- 
pling Jij by J[j = Jij + SKij . Here ATy is again a Gaus- 
sian distributed random number with variance one and 
the parameter S measures the strength of the perturba- 
tion. The comparison of the energy balance Ai^dofoct for 
turning over a connected spin cluster of linear extent L 
with the change of the ground state energy AE'random 
induced by the random variation of the couplings yields 
an estimate for the length scale beyond which the origi- 
nal ground state is unstable with respect to the pertur- 
bation. Ai?dcfoct is simply the defect energy, which is 
proportional to JL^ (see section ||). The contribution 
to AE'random comiug from the L"^^ interface spins of the 
cluster {ds being the fractal dimension of the interface) 
is proportional to 5L'^^/'^. Thus for L > L*{5) with 



L*{5)^{JI5) 



-i/C 



with C = ds/2-y (25) 



we have AA^randomCi) > Ai?dofGct(-^) and flipping of clus- 
ters is favored by the perturbation. Thus the ground 
state configurations in the original (denoted by Si) and 
the perturbed sample (denoted by Sl{S)) become uncor- 
related for distances larger than L* . 

This statement can be quantified by studying the over- 
lap correlation function 



Csir) 



1 ^ 

— S, S^+r Sl{S) 5-+^ ((5) 



(26) 



According to the above mentioned argument one expects 
in the limit A^ — > oo a scaling form 



Cs{r) r^£{rS^/'^). 



(27) 



In figure ^ we show the result of our calculation of the 
overlap correlation function Cs{r). We fixed the system 
size to L = 50, for which reason one has to neglect the 
data points for r > L/4 (note the upwards bending due 
to the periodic boundary conditions). For the rest of the 
data we obtain the best data collapse for 

1/C = 1.05 ±0.05 i.e. C = 0.95 ± 0.05 (28) 

which j-agrees well with the estimate from Bray and 
MooreE3 obtained in a different way and by considering 
smaller system sizes. With the value for y we reported 
in section || the fractal dimension of the interface of an 
excitation is given by ds = 1.34 ± 0.10. 

In passing we mention that the dependency of C(r) 
on distance r is neither exponential nor algebraic: it can 
nicely be fitted with a stretched exponential 



C{r) « exp(-r''/6) + exp(-(L - r)76) 



(29) 
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FIG. 6. Scaling plot of the overlap correlation function 
Cs{r) versus r/L* with L* = (J"^/"-. The best data collapse 
(for data confined to r < L/4) is obtained for 1/^ = 1.05. 
The system size is L = 50 and the data are averaged over 
400 samples. These were obtained by creating > 80 reference 
instances and creating 5 random perturbations of strength 5 
for each. 



with fit parameters a and b. For instance 5 = 0.1 for 
L = 50 yields a = 0.8 and b = 0.75. Since a and b 
seem to depend on the perturbation strength S we do 
not expect the form ( ^ ) to be universal. 

Defining fL((5) — J2r=Q^s{r) one expects from ( p7| ) 

US)^i{LS'^<) (30) 
and in a more direct way for the ground state overlapEl 

QLi6) = \EliS,siiS)\ 



Ql{5) = Q{L5^'^) 



(31) 



Note that — Q\{5). We show a finite size scahng 

plot for Ql{S) in fig. 0, from which we estimate 1/C = 
1.2 ± 0.1. The quality of the data collapse is good (cf. 
fig. 2 of refE3). 



VI. SUMMARY 

With the help of an improved branch & cut algorithm 
we were able to reinvestigate the critical behavior of the 
two-dimensional Ising spin glass model with a continuous 
bond distribution with much better accuracy. We found 
that the stiffness exponent is given hy y = —0.281 ±0.002 
implying a correlation length exponent olv — 3.56±0.02, 
which agrees well with our independent estimate v — 
3.6 ±0.2 from Monte Carlo simulations. For the latter we 
introduced a pseudo-ferromagnetic order parameter with 
the help of exactly known ground states and analyzed its 
finite size scaling behavior at non-zero temperatures. 
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FIG. 7. Scaling plot of the ground state overlap Ql{S). 
The best data collapse is obtained with 1/C = 1-2. 



We hope that our calculation settles the controversy re- 
garding the thermal exponejUt v initiated by the cluster 
Monte Carlo study of LiangH and tfae numerical transfer 
matrix study by Kawashima et alB: Their values for i/ 
are substantially smaller than ours indicatina-a violation 
of the scaling prediction by Bray and MooreO. Our re- 
sults for y and are clearly compatible with this scaling 
prediction v = —l/y. 

Furthermore we determined exact ground states for 
systems within an external field and from a finite size 
scaling analysis of the magnetization we obtained an 
independent estimate for the magnetic exponeijtt, 6 ~ 
1.48 ± 0.01. This confirms an earlier observationtS that 
there seems to be a disagreement between the scaling 
theoryO predicting b = \ — y and the numerical values 
obtained so far. In particular this discrepancy does not 
fade away for larger system sizes, which we were able to 
study here. Therefore our conclusion is that there must 
be a deeper reason for this disagreement than some finite 
size effect which might disappear if one only considers 
large enough system sizes. 

Moreover, we calculated the overlap correlation func- 
tion by perturbing the bonds slightly in a random man- 
ner. We found a chaos exponent C — 0.95 ±0.05 in agree- 
ment with earlier estimates from the analysis of smaller 
system sizes. 

Finally a few words concerning future perspectives: 
First we would like to point out that in principle it is 
possible to improve the system sizes and quality of statis- 
tics even further with the algorithm we have at hand, 
provided we could simply run it on a powerful paral- 
lel machine. However, our algorithm relies heavily on 
a commercial linear problem solver for which we do not 
have a license to run it on hundreds of processors of a e.g. 
Paragon XP/SIO. On such a machine we could possibly 
obtain an acceptable quality of statistics for L = 100, for 
which we can presently do only a few samples in reason- 
able time on individual workstations. 

As has been mentioned in the introduction, recently a 
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finite temperature phase transition iii the site disordered 
Ising spin glass has been reportedo. Since the critical 
temperature is pretty small, though, Monte Carlo studies 
might be hampered by equilibration problems. Therefore 
this result could be put on a much firmer base, if the 
stiffness exponent y would indeed turn out to be positive 
in this particular two-dimensional model and so signaling 
the stability of the spin glass ordered phase for small, 
non-vanishing temperatures. We intend to answer this 
question with our algorithm soon. 

Furthermore an obvious and highly rewarding step 
would be to perform the same study in three dimensions. 
To calculate ground states for the three-dimensional Ising 
spin glass model is an NP-complete problem and the 
two-dimensional problem we have studied here is NP- 
complete, too (note we have a continuous bond distri- 
bution, periodic boundary conditions and an external 
field). However, although both questions belong to the 
same class of hard combinatorial problems the three- 
dimensional Ising spin glass is much harder, which means 
that the operation count will be much higher: either the 
power of the L, the system size, or the coefficient in the 
exponent will be larger for three dimensions than for two 
dimension. Nevertheless we are currently undertaking ef- 
forts in this direction, our progress in this matter will be 
reported elsewhere. 

ACKNOWLEDGMENTS 

We should like to thank A. Bray, C. De Simone, G. 
Reinelt, G. Rinaldi, D. Stauffer, and A. P. Young for 
various stimulating discussions and valuable hints and 
remarks. The Monte Carlo simulations have been done 
on the Paragon XP/SIO of the ZAM at KFA Jiilich. 



* W. L. McMillan, Phys. Rev. B 29, 4026 (1984); A. J. Bray 
and M. A. Moore , J. Phys. C 17, L463 (1984); ibtd. L613. 

® N. Kawashima, N. Hatano and M. Suzuki, J. Phys. A 25, 
4985 (1992). 

" I. Bieche, R. Maynard, R. Rammal and J. P. Uhry, J. Phys. 
A 13, 2553 (1980). 

^ R. M. Karp, p. 85 in Complexity of Computer Computa- 
tions ed. R. E. Miller and J. W. Thatcher (Plenum Press, 
New York 1972). 

^ F. Barahona, J. Phys. A 15, 3241 (1982); Phys. Rev. B 49, 
12864 

^ M. Grotschel, M. Jiinger and G. Reinelt, in: Heidelberg 
Colloqium on Glassy dynamics and Optimization , ed. L. 
van Hemmen and I. Morgenstern (Springer- Verlag, Heidel- 
berg 1985); 

" C. De Simone, M. Diehl, M. Jiinger, P. Mutzel, G. Reineh 
and G. Rinaldi, J. Stat. Phys. 80, 487 (1995); ibtd. 84, 
(1996). 

^ L. Saul and M. Kardar, Nuc. Phys. B 432 641 (1994). 

^ N. Kawashima and M. Suzuki, J. Phys. A 25, 1055 (1992). 

^ A. J. Bray and M. A.Moore in: Heidelberg Colloqium on 
Glassy dynamics and Optimization, ed. L. van Hemmen 
and I. Morgenstern (Springer- Verlag, Heidelberg 1985); in: 
Advances on Phase Transitions and Disorder Phenomena, 
ed. G. Busiello and L. Deces, (World Scientific, Singapore 
1986). 

* F. Barahona, M. Grotschel, M. Jiinger and G. Reinelt, Op- 
erations Research 3, 493 (1988). 

^ Bray, A. J., Comments Cond. Mat. Phys., 14 21 (1988). 
J. Kisker, L. Santen, M. Schreckenberg and H. Rieger, 
Phys. Rev. B 52 xxx (1996). 

M. Jiinger and G. Rinaldi, Technical-Report (in prepara- 
tion). 

A. J.Bray and M. A.Moore, Phys. Rev. Lett. 58, 57 (1987). 
'■^ Y. Imry and S.-k. Ma, Phys. Rev. Lett. 35, 1399 (1975). 



^ For a review see H. Rieger Monte Carlo studies of Ising 
spin glasses and random field systems, Ann. Rev. Comp. 
Phys. II, p. 296, ed. D. Stauffer (World Scientific, 1995) 
and references therein. 

^ For site-disorder instead of bond disorder there are indi- 
cations of a finite temperature phase transition in 2d: T. 
Shirakura and F. Matsubara, J. Phys. Soc. Jap. 64 2338 
(1995). 

3 C. Dekker, A. F. M. Arts and H. W. de Wijn, Phys. Rev. 

B 38, 11512 (1988). 
* R. N. Bhatt and A. P. Young, Phys. Rev. B 37 5606 (1988). 
^ S. Liang, Phys. Rev. Lett. 69, 2145 (1992). 
® R. R. P. Singh and S. Chakravarty, Phys. Rev. Lett. 57, 

245 (1986). 

^ H. F. Cheung and W. L. McMillan, J. Phys. C 16, 7033 
(1983); D. Huse and I. Morgenstern, Phys. Rev. B 32, 3032 
(1985). 



8 



